# construction network

# Rscript step1.R --file "../test_data/datExpr.csv" --output "../1_construct_network"

library(optparse)


option_list = list(
  make_option("--exprFile", type="character", default=NULL,
              help="input file name"),
  make_option("--output", type="character", default="out.txt",
              help="output file name [default= %default]"),
  make_option("--TOMType", type="character", default="unsigned",
              help="TOMType [default= %default]"),
  make_option("--minModuleSize", type="integer", default=30,
              help="minModuleSize [default= %default]"),
  make_option("--maxPower", type="integer", default=20,
              help="maxPower [default= %default]"),
  make_option("--mergeCutHeight", type="double", default=0.25,
              help="mergeCutHeight [default= %default]")
)

# 解析命令行参数
opt_parser = OptionParser(option_list=option_list, add_help_option=TRUE)
opts = parse_args(opt_parser)

datExpr = read.csv(opts$exprFile, stringsAsFactors = FALSE, row.names = 1);
nGenes = ncol(datExpr);
nSamples = nrow(datExpr);

library(WGCNA)


# auto select soft-thresholding power
# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=ifelse(opts$maxPower > 20, opts$maxPower, 20), by=2))
# Call the network topology analysis function
allowWGCNAThreads()
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

png(
  filename = paste(opts$output, "chose_power.png", sep="/"),
  width = 9,
  height = 5,
  units = "in",
  res = 300
)

par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

dev.off()


soft_power = sft$powerEstimate

# 评估网络的连通性
png(
  filename = paste(opts$output, "check_topology.png", sep="/"),
  width = 9,
  height = 5,
  units = "in",
  res = 300
)

adjacency = abs(cor(datExpr, use="p"))^soft_power
k = as.vector(apply(adjacency,2,sum,na.rm=T))

par(mfrow=c(1,2))
cex1=0.9

hist(k)
scaleFreePlot(k)
dev.off()


net = blockwiseModules(
  datExpr,
  power = soft_power,
  TOMType = opts$TOMType,
  minModuleSize = opts$minModuleSize,
  reassignThreshold = 0,
  mergeCutHeight = opts$mergeCutHeight,
  numericLabels = TRUE,
  pamRespectsDendro = FALSE,
  saveTOMs = FALSE,
  verbose = 3
)

TOM = TOMsimilarityFromExpr(datExpr, power = soft_power);

# save net data
save(net, file=paste(opts$output, "network.RData", sep="/"))
save(TOM, file=paste(opts$output, "TOM.RData", sep="/"))

# sample cluster
sampleTree = hclust(dist(datExpr), method = "average");
sample_num = nrow(datExpr)

png(
  filename = paste(opts$output, "sample_cluster.png", sep="/"),
  width = 13,
  height = 9,
  units = "in",
  res = 300
)

plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)

dev.off()
